www.gusucode.com > 溷沌分析工具箱 - Chaotic Systems Toolbox > 混沌分析工具箱 - Chaotic Systems Toolbox\noisergeo.m
function [xr,Yr]=noisergeo(x,dim,tau,r,q,p,theta)
%Syntax: [xr,Yr]=noisergeo(x,dim,tau,r,q,p,theta)
%________________________________________________
%
% Noise reduction by Local Geometric Projection.
%
% xr is the vector/matrix with the cleaned time series.
% Yr is the phase space of the last cleaned xr.
% x is the time series.
% dim is the embedding dimension.
% tau is the time delay.
% r can be either
% real defining the neighborhood range.
% integer defining the number of nearest neighbors.
% q can take one of the following values
% 'wAV' for the weighted average.
% [an integer from 0 to dim-1] for the local geometric projection.
% 'mod' for the adaptive selection of the local neighborhood dimensions.
% p defines the norm.
% theta is the correction paprameter [0,1]
% 1: full correction.
% 0: no correction.
%
%
% References:
%
% Kantz H, Schreiber T, Hoffmann I, Buzug T, Pfister G, Flepp L G, Simonet
% J, Badii R, Brun E (1993): Nonlinear noise reduction: A case study on
% experimental data. Physical Review E 48: 1529-1538
%
% Leontitsis A., Bountis T., Pange J. (2004): An adaptive way for improving
% noise reduction using local geometric projection. CHAOS 14(6): 106-110
%
%
% Alexandros Leontitsis
% Department of Education
% University of Ioannina
% Ioannina
% Greece
%
% University e-mail: me00743@cc.uoi.gr
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% 14 Jul 2001
if nargin<1 | isempty(x)==1
error('You should provide a time series.');
else
% x must be a vector
if min(size(x))>1
error('Invalid time series.');
end
x=x(:);
% n is the time series length
n=length(x);
end
if nargin<2 | isempty(dim)==1
dim=2;
else
% dim must be either a scalar or a vector
if min(size(dim))>1
error('dim must be a scalar or a vector.');
end
% dim must be an integer
if round(dim)-dim~=0
error('dim must be an integer.');
end
% dim values must be above 1
if any(dim<1)==1
error('dim values must be above 1');
end
end
if nargin<3 | isempty(tau)==1
tau=1;
else
% tau must be either a scalar or a vector
if min(size(tau))>1
error('tau must be a scalar or a vector.');
end
% tau must be an integer
if round(tau)-tau~=0
error('tau must be an integer.');
end
% tau values must be above 1
if any(tau<1)==1
error('tau values must be above 1');
end
end
if nargin<4 | isempty(r)==1
r=dim+1;
else
% r must be either a scalar or a vector
if min(size(r))>1
error('r must be a scalar or a vector.');
end
% r values must be above 0
if any(r<=0)==1
error('r values must be greater than 0');
end
end
if nargin<5 | isempty(q)==1
q=1;
end
if nargin<6 | isempty(p)==1
p=2;
else
% p must be either a scalar or a vector
if min(size(p))>1
error('p must be a scalar or a vector.');
end
end
if nargin<7 | isempty(theta)==1
theta=1;
else
% theta must be either a scalar or a vector
if min(size(theta))>1
error('theta must be a scalar or a vector.');
end
% theta must be above 0 and bellow 1
if theta<0 | theta>1
error('theta must be above 0 and bellow 1.')
end
end
% Only one of dim, tau, r, p or theta should be vector
l=[length(dim),length(tau),length(r),length(p),length(theta)];
if length(find(l>1))>1
error('Only one of dim, tau, r, p, or theta should be vector.');
end
% Make the phase-space
[Y,T]=phasespace(x,dim,tau);
m=max(l);
dim=ones(1,m).*dim;
tau=ones(1,m).*tau;
r=ones(1,m).*r;
p=ones(1,m).*p;
theta=ones(1,m).*theta;
for i=1:m
% Initialize Yr
Yr=zeros(T,dim(i));
% For every phase-space point
for j=1:T
% Locate the j-th point
y=Y(j,:);
% Check neighborhood or neighbors
if mod(r(i),floor(r(i)))==0
lock=Knearest(y,Y,r(end)+1,p(i));
else
lock=radnearest(y,Y,T,r(i),p(i));
lock(find(j==lock))=[];
end
if isempty(lock)==1
Yr(j,:)=y;
elseif q=='wAV'
% The calculations for the weighted average
Ynearest=Y(lock,:);
w=[];
for k=1:length(lock)
w(k)=norm(y-Ynearest(k,:))/2/r(i);
w(k)=exp(-w(k)^2);
Ynearest(k,:)=Ynearest(k,:)*w(k);
end
Yr(j,:)=sum(Ynearest)/sum(w);
else
% All the neighboring points have equal weight
if length(lock)<dim(i)
Yr(j,:)=mean(Y(lock,:));
else
% Make the matrix of the local neighborhood
Ynearest=Y(lock,:);
%Ynearest=quality3(Ynearest,2);
%Ynearest=quality2(Ynearest,2,.99,r,y);
Ynearest=quality(Ynearest,2,.7);
% Neigborhood average
ybar=mean(Ynearest);
% Subtract the average from the neighborhood poionts
Ynearest=Ynearest-ones(size(Ynearest,1),1)*ybar;
% SVD on the local neighborhood
[u,s,v]=svd(Ynearest,0);
% Either use the Adaptive selection ...
if q=='mod'
s=diag(s);
d=s(1:dim(i)-1)./(s(2:dim(i))+eps);
q=find(d==max(d));
Yr(j,:)=y-(y-ybar)*(v(:,q+1:dim(i))*v(:,q+1:dim(i))');
% ... or use a predifined q
else
Yr(j,:)=y-(y-ybar)*(v(:,dim(i)-q:dim(i))*v(:,dim(i)-q:dim(i))');
end
end
end
Yr(j,:)=theta(i)*Yr(j,:)+(1-theta(i))*Y(j,:);
end
% Initialize xr
xr(:,i)=zeros(n,1);
times(:,i)=zeros(n,1);
% Calculate the reconstructed time series
for j=1:dim
xr(1+(j-1)*tau:T+(j-1)*tau,i)=xr(1+(j-1)*tau:T+(j-1)*tau,i)+Yr(:,j);
times(1+(j-1)*tau:T+(j-1)*tau,i)=times(1+(j-1)*tau:T+(j-1)*tau,i)+1;
end
xr(:,i)=xr(:,i)./times(:,i);
xr(:,i)=xr(:,i)-mean(xr(:,i)-x);
end